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a  pattern  of  maxima  and  minima  on  the  line  normal  to  the  path  passing  through 
the  center  of  the  disturbance.  As  expected,  the  focusing/diffraction  effects 
diminish  when  the  transverse  dimension  of  the  disturbance  exceeds  the  width  of 
the  first  Fresnel  zone — typically,  several  megameters. 

The  analysis  models  widespread  inhomogeneities,  such  as  within  the  polar 
cap  or  at  the  day/night  terminator,  as  semiinfinite  regions  separated  by  dif¬ 
fuse  boundaries;  it  then  derives  full-wave  analytic  expressions  for  the  re¬ 
flection  of  the  TEM  mode.  Mode  reflection  is  found  to  significantly  affect  an 
ELF  signal  in  two  actual  situations:  first,  when  receivers  are  on  great  circle 
paths  that  are  nearly  tangential  to  the  disturbed  polar  cap — in  which  case 
shadow  zones  and  interference  patterns  can  occur;  and  second,  when  signals  are 
incident  on  the  day/night  terminator  (from  the  day  side)  at  angles  exceeding 
about  75  deg — in  which  case  the  signals  are  affected  by  a  phenomenon  analogous 
to  total  internal  reflection.  Reflection  is  found  to  be  unimportant  if  the 
boundary  thickness  exceeds  about  one-sixth  of  a  wavelength. 
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PREFACE 


This  report  concludes  Pacific-Sierra  Research  Corporation's  (PSR's) 
analysis  of  the  effects  of  nonstratif ied  ionospheric  disturbances  on 
ELF  propagation.  It  extends  the  treatments  presented  in  two  previous 
documents:  ELF  Propagation  in  a  Non-Stratified  Earth  Ionosphere 
Waveguide ,  PSR  Report  806,  April  1978;  and  An  Integral  Equation 
Approach  to  ELF  Propagation  in  a  Non-Stratified  Earth- Ionosphere 
Waveguide 3  PSR  Report  904,  February  1979.  The  work  was  sponsored  by 
the  Office  of  Naval  Research,  Arlington,  Virginia. 
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SUMMARY  AND  CONCLUSIONS 

This  report  analyzes  the  propagation  of  the  TEM-ELF  waveguide 
mode  when  the  ionosphere  is  not  stratified.  It  treats  strong  local¬ 
ized  ionospheric  disturbances  by  recasting  the  lateral  wave  equation 
as  a  two-dimensional  integral  equation,  and  applies  a  specially  de¬ 
veloped  algorithm  to  obtain  numerical  solutions.  The  quasi- full-wave 
results  show  that  a  localized  ionospheric  disturbance  behaves  like  a 
converging  cylindrical  lens  filling  a  narrow  aperture.  Lateral  dif¬ 
fraction  and  focusing,  ignored  in  treatments  that  do  not  fully 
account  for  transverse  ionospheric  structure,  cause  the  ELF  signal 
to  exhibit  a  pattern  of  maxima  and  minima  on  the  line  normal  to  the 
path  passing  through  the  center  of  the  disturbance.  As  expected, 
the  focusing/diffraction  effects  diminish  when  the  transverse  dimen¬ 
sion  of  the  disturbance  exceeds  the  width  of  the  first  Fresnel  zone — 
typically,  several  megameters. 

The  analysis  models  widespread  inhomogeneities,  such  as  within 
the  polar  cap  or  at  the  day/night  terminator,  as  semiinfinite  regions 
separated  by  diffuse  boundaries;  it  then  derives  full-wave  analytic 
expressions  for  the  reflection  of  the  TEM  mode.  Mode  reflection  is 
found  to  significantly  affect  an  ELF  signal  in  two  actual  situations: 
first,  when  receivers  are  on  great  circle  paths  that  are  nearly  tan¬ 
gential  to  the  disturbed  polar  cap — in  which  case  shadow  zones  and 
interference  patterns  can  occur;  and  second,  when  signals  are  incident 
on  the  day/night  terminator  (from  the  day  side)  at  angles  exceeding 
about  75  deg — in  which  case  the  signals  are  affected  by  a  phenomenon 
analogous  to  total  internal  reflection.  Reflection  is  found  to  be 
unimportant  if  the  boundary  thickness  exceeds  about  one-sixth  of  a 
wavelength. 


» 


PRECEDING  PAGE  BUNK  -  NOT 


FLAMED 


-vii- 


CONTENTS 


PREFACE  .  ill 

SUMMARY  AND  CONCLUSIONS  .  v 

FIGURES  .  ix 


Section 

I.  INTRODUCTION  .  1 

II.  PROPAGATION  IN  LATERALLY  NONUNIFORM  EARTH- 

IONOSPHERE  WAVEGUIDE  .  3 

III.  MODELS  OF  WAVEGUIDE  DISTURBANCES  .  7 

Laterally  Bounded  Disturbance  .  7 

Widespread  Disturbance  .  7 

IV.  REFLECTION  COEFFICIENT  .  10 

V.  NUMERICAL  RESULTS  .  14 

Semiinfinite  Disturbed  Region  .  14 

Bounded  Disturbance  .  23 

APPENDIX .  33 


REFERENCES 


37 


\\  PBECSDIMG  pace  blank  -  N€  r 

_  mil  —iiiiMiir - -•>"  'Y 


► 


-IX- 


FIGURES 


1.  Schematic  of  propagation  for  laterally  bounded 

ionospheric  disturbance  . . .  6 

2.  Schematic  of  reflection  from  disturbed  region  with 

diffuse,  stratified  boundary  .  8 

3.  Reflection  from  strong  SPE  (F  =  45  Hz)  .  16 

4.  Reflection  from  strong  SPE  (F  =  75  Hz)  .  17 

5.  Reflection  from  moderate  SPE  (F  =  45  Hz)  .  18 

6.  Reflection  from  weak  SPE  (F  =  45  Hz)  .  19 

7.  Reflection  from  day/night  terminator:  mode  incident 

from  nighttime  side  (F  =  45  Hz)  .  20 

8.  Reflection  from  day/night  terminator:  mode  incident 

from  daytime  side  (F  =  45  Hz)  .  21 

9.  Amplitude  of  W  versus  off-path  distance  y  for  x  = 

10  Mm,  x  =  5  Mm,  Ax  =  0.5  Mm,  and  strong  (A  =  1) 
disturbance  with  various  half-widths  .  25 

10.  Amplitudeof  W  versus  off-path  distance  y  for  x  = 

10  Mm,  x  =  5  Mm,  Ay  =  Ax  =  0.5  Mm,  and  various 

disturbance  strengths  .  27 

11.  On-path  amplitude  of_W  versus  disturbance  half-width 

Ay  for  x  =  10  Mm,  x  =  5  Mm,  Ax  =  0.5  Mm,  y  *  0, 

and  various  disturbance  strengths  .  28 

12.  Amplitude  of  transverse  shape  factor  Q  versus  dis¬ 

turbance  half-width  Ay  for  disturbance  centered  on 

direct  propagation  path  (y  *  0)  .  30 

13.  Amplitude  of  W  versus  distance  along  direct  propaga¬ 

tion  path  for  strong  disturbance  (A  =  1,  y  =  0, 

Ay  =  1  Mm)  .  31 

A.l.  Partition  for  solution  within  disturbed  region  .  34 

A. 2.  Partition  for  solution  outside  disturbance  .  34 


I 


t 


» 


-1- 


I.  INTRODUCTION 


This  report  treats  the  propagation  of  the  transverse  electro¬ 
magnetic  (TEM)  extremely- low-frequency  (ELF)  waveguide  mode  under 
conditions  in  which  the  ionosphere  is  not  stratified.  We  consider 
both  widespread  and  laterally  bounded  inhomogeneities  of  arbitrary 
strength,  and  impose  no  restrictions  on  the  direction  of  propagation 
relative  to  ionospheric  gradients.  The  report  therefore  supplements 
Field  and  Joiner's  [1979]  integral-equation  analysis  of  ELF  propaga¬ 
tion,  which  yielded  numerical  results  for  (1)  weak,  localized  iono¬ 
spheric  disturbances  remote  from  both  transmitter  and  receiver,  and 
(2)  ionospheric  disturbances  of  arbitrary  strength  and  extent,  azi- 
muthally  symmetric  about  the  transmitter.  Findings  accounted,  in  the 
first  case,  for  gradients  transverse  to  the  propagation  path,  and  in 
the  second  case,  for  gradients  parallel  to  the  direction  of  propaga¬ 
tion.  The  results  also  showed  the  relationship  between  the  fraction 
of  the  first  Fresnel  zone  filled  by  an  ionospheric  disturbance  and 
the  resulting  propagation  anomaly. 

To  analyze  the  propagation  effects  of  widespread  ionospheric  inho¬ 
mogeneities,  the  present  analysis  develops  a  computational  model  that 
assumes  two  laterally  uniform — but  vertically  nonuniform — semiinfinite 
regions.  The  model  is  adequate  to  describe  reflections  from  such  wide 
areas  as  the  polar  cap  during  a  solar  proton  event  (SPE),  the  zone 
affected  by  a  solar  X-ray  flare,  and  the  day  and  night  portions  of 
the  earth,  and  would  also  apply  to  the  disturbance  created  by  a 
nuclear  burst  at  altitudes  above  a  few  hundred  kilometers.  Bounding 
such  areas  are  transition  regions  in  which  the  waveguide  varies  over 
a  few  hundred  to  perhaps  two  thousand  kilometers.  Since  those  dis¬ 
tances  compare  to  A/2tt  (where  X  is  wavelength),  the  simplifications 
associated  with  either  sharp  (Fresnel  coefficients)  or  gentle  (Wenzel- 
Kramers-Brillouln  [WKB]  calculations)  boundaries  cannot  be  used  at 
ELF.  We  therefore  use  a  full-wave  theory  that  completely  accounts 
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for  the  lateral  diffuseness  of  the  boundaries  but  still  yields 

* 

closed-form  solutions  for  the  reflection  coefficients. 

Closed-form  solutions  are  not  possible  for  strong  ionospheric 
disturbances  that  are  Loo  localized  for  approximation  as  semiinfinite 
regions.  To  treat  the  localized  situation,  we  recast  the  lateral 
wave  equation  as  a  two-dimensional  integral  equation,  then  develop  a 
special  algorithm  to  obtain  numerical  solutions.  The  full-wave  re¬ 
sults  account  for  the  lateral  focusing  and  diffraction  omitted  from 
earlier,  less-exact  formulations. 

Section  II  reviews  the  two-dimensional  wave  equations  that  ap¬ 
proximately  describe  ELF  propagation  in  a  laterally  nonuniform  earth- 
ionosphere  waveguide.  Section  III  presents  models  for  the  waveguide 
propagation  constant  in  ambient  and  disturbed  regions.  Section  IV 
adapts  the  Epstein  theory  of  propagation  in  vertically  inhomogeneous 
ionospheres  so  as  to  yield  the  lateral  reflection  coefficient  of  the 
TEM-ELF  waveguide  mode.  Section  V  presents  numerical  results  for 
several  model  disturbances.  The  Appendix  describes  the  algorithm 
used  to  solve  the  two-dimensional  integral-wave  equation. 


The  theory  applies  only  to  ELF,  which  is  characterized  by  a 
single  propagating  waveguide  mode.  At  higher  frequencies,  any  bound¬ 
ary  abrupt  enough  to  require  the  full-wave  theory  would  cause  coupling 
among  higher  order  modes,  not  accounted  for  here. 


II.  PROPAGATION  IN  LATERALLY  NONUNIFORM  EARTH-IONOSPHERE  WAVEGUIDE 


The  wave  equation  can  be  solved  by  separation  of  variables  under 

* 

conditions  in  which  the  ionosphere  is  laterally  uniform.  The  re¬ 
sult,  well  known  from  Galejs  [1972],  for  example,  is 


Eq  =  AAqFq(z)Eq(x,  y)  V/m  , 


(1) 


where  is  the  vertical  component  of  the  electric  field;  A  is  a  con¬ 
stant  involving  dipole  moment,  wave  frequency,  and  ground  conductivity 
Aq  is  the  excitation  factor  representing  the  efficiency  with  which  the 
TEM  waveguide  mode  is  launched;  F^(z)  is  the  height-gain  function  for 
the  vertical  dependence  of  the  field,  normalized  to  unity  at  z  =  0; 
x,  y,  and  z  are  Cartesian  coordinates;  and  the  subscript  0  denotes 
undisturbed,  laterally  homogeneous  conditions. 

The  lateral  dependence  of  the  field  is  governed  by  the  function 
EQ,  which  satisfies  the  two-dimensional  wave  equation 


(vt  +  ,Ao)V°  •  <2) 

2 

where  is  the  two-dimensional  Laplacian,  k  is  the  free-space  wave 
number,  and  Sq  is  a  waveguide  propagation  constant  determined  by  im¬ 
posing  boundary  conditions  on  Fq  at  the  ground  and  in  the  ionosphere. 
The  literature  supplies  full-wave  methods  for  calculating  Sq  and  A^ 
for  virtually  any  ionospheric  height  profile,  as  well  as  numerical  re¬ 
sults  for  many  models  of  ambient  and  disturbed  ionospheres  [Budden, 
1961;  Field,  1970;  Wait,  1970;  Galejs,  1972;  Pappert  and  Moler,  1974; 
Greifinger  and  Griefinger,  1978].  Once  Sq  is  determined,  Eq.  (2)  is 
easily  solved  for  E^. 


To  avoid  complexities  unrelated  to  lateral  ionospheric  gradients 
our  calculations  ignore  the  earth's  curvature — an  approximation  well 
justified  for  ELF  propagation  parameters  other  than  spreading  loss. 
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To  characterize  in  physical  terms,  we  note  that  at  large  lateral 
distances  D  from  the  source, 

EQ  oc  exp  (-ikS^D)  .  (3) 

Thus,  ReS^  is  the  ratio  of  the  speed  of  light  to  the  phase  velocity  of 
the  waveguide  mode,  and  ImS^  is  proportional  to  the  attenuation  coef¬ 
ficient  a,  such  that 

a  =  -8.6kImSQ  dB/Mm  ,  (4) 

where  k  is  in  Mm  . 

In  the  presence  of  lateral  ionospheric  gradients,  a  rigorous 
separation  of  variables  is  impossible.  Hence,  fields  cannot  be  ex¬ 
pressed  as  products  of  vertical  and  radial  functions,  as  they  are  in 
Eq .  (2).  But  when  ionospheric  irregularities  are  large  enough  to  sig¬ 
nificantly  affect  ELF  propagation,  scales  for  lateral  variations  of 
the  ionospheric  refractive  index  tend  to  be  at  least  ten  times  longer 
than  those  for  vertical  variations.  It  can  therefore  be  argued  that 
the  waveguide  propagation  parameters  depend  mainly  on  the  local  iono¬ 
sphere,  and  that  separation  of  the  fields  into  vertical  and  lateral 
functions  is  an  approximately  correct  procedure. 

Using  the  above  argument.  Field  and  Joiner  [1979]  combined 
eikonal  and  full-wave  methods  to  calculate  the  fields  in  a  laterally 
nonuniform  waveguide.  Their  main  approximation  was  to  separate  the 
ground-level  fields  into  vertical  and  lateral  functions,  each  of  which 
was  then  computed  on  a  full-wave  basis.  The  ground-level  field  [note 
that  F(z  =  0)  =  1]  was  assumed  to  be 

E  «  AA(x,  y)E(x,  y)  ,  (5) 

/v 

where  the  lateral  dependence  resides  mainly  in  E,  satisfying 

[v2  +  k2S2(x,  y)]  E  =  0  .  (6) 
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The  full-wave  methods  for  obtaining  Sq  and  Aq  (given  in  the  ref¬ 
erences  cited  above)  can  also  be  used  to  obtain  A(x,  y)  in  Eq.  (5)  and 
S(x,  y)  in  Eq.  (6).  Such  calculations  must  be  carried  out  for  many 
locations  (x,  y),  each  governed  by  a  local  ionospheric  height  profile, 
to  well  represent  the  lateral  dependence  of  S  and  A.  A  single  calcula¬ 
tion  of  Sq  and  Aq,  on  the  other  hand,  suffices  to  represent  all  loca¬ 
tions  under  laterally  uniform  conditions. 

To  calculate  reflections  from  widespread  disturbances  that  can  be 
approximated  as  semiinfinite  regions  with  diffuse  boundaries,  we 
solve  Eq.  (6)  directly  (see  Sec.  IV).  To  calculate  the  propagation 
anomaly  due  to  a  laterally  bounded  disturbance,  however,  it  is  con¬ 
venient  to  recast  the  problem  as  an  integral  equation.  A  laterally 

nonuniform  ionospheric  disturbance  can  be  characterized  by  the  dif- 

2  2  /*. 
ference  S  (x,  y)  -  Sg.  The  undisturbed  wave  function  Eq  is  governed 

by  Sq  and  can  be  assumed  known.  Subtracting  Eq.  (2)  from  Eq.  (4)  gives 
(V2  +  k2S*)(ff  -  Eg)  =  -k2(S2  -  S2)E  .  (7) 


t 


Then,  using  the  Green's  function 

G  =  -i7TH^2)  (kSQr2)  ,  (8) 

where  H  is  the  Hankel  function,  and  applying  the  two-dimensional  Green's 
theorem,  the  integral  equation  for  Ef  is  obtained: 

E(x,  y)  =  Eg(x,  y) 

-  nh  ff  *’  dy'  [s2(x’’  ,')  '  so] 

— OO 

•  Hj2)(kS0r2)E(x',  y')  .  (9) 

Equation  (9)  is  formally  identical  to  an  integral  equation  of  Wait 
[1964].  The  diagram  of  the  propagation  in  Fig.  1  defines  the  geo¬ 
metric  terms  in  the  equations  above. 
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Integration  point 


Receiver 

(x,y) 


Fig.  1- -Schematic  of  propagation  for  laterally  bounded 
ionospheric  disturbance 


It  is  also  convenient  to  define  a  relative  propagation  function  W, 
denoting  the  fractional  amount  by  which  the  disturbed  lateral  wave 
function  ff  differs  from  the  undisturbed  function  E^:  specifically, 

lf(x,  y)  J  W(x,  y ) e’q ( x ,  v)  .  (10) 

Inserting  Eq .  (10)  into  Eq .  (9)  and  expressing  as  for  a  horizontal 
electric  dipole  transmitter  [Field  and  Joiner,  1979]  gives 


00 


(kSQr2)H1 

H<2)(kSQr) 


W(x\  y')  , 


(11) 


which  is  the  most  general  fora  of  the  integral  equation  for  the  rela¬ 
tive  propagation  function. 


I II.  MODELS  OF  WAVEGUIDE  DISTURB AN CES 


The  analytic  models  given  below  pertain  to  two  types  of  wave¬ 
guide  disturbances — laterally  bounded  and  widespread,  both  of  arbi¬ 
trary  strength — that  were  not  analyzed  by  Field  and  Joiner  [1979]. 

The  numerical  results  reported  in  Sec.  V  were  obtained  with  the 
present  models. 

LATERALLY  BOUNDED  DISTURBANCE 

A  laterally  bounded  disturbance  was  schematically  illustrated 
in  Fig.  1.  We  represent  the  disturbed  region  by 

[.?  -  «»  [-  -p  [-  .  u»> 

2 

where  denotes  the  value  of  the  waveguide  propagation  constant  at 
x  =  x,  y  =  0.  Note  that  Eq.  (12)  implies  no  limits  on  the  strength 
Sj  -  Sgj  and  that  the  disturbance  is  assumed  to  center  on  the  x-axis, 
with  off-path  effects  accounted  for  by  adjusting  the  transverse  loca- 

2  2  2  2  — 

tion  y  of  the  receiver.  Adjusting  -  Sg,  (Ax)  ,  (Ay)  ,  and  x  hence 
allows  the  equation  to  represent  disturbances  with  different  strengths, 
lateral  gradients,  and  longitudinal  positions.  The  disturbance  effec¬ 
tively  vanishes  when  |x  -  xj  or  y  becomes  several  times  greater  than 
Ax  or  Ay. 

WIDESPREAD  DISTURBANCE 

The  general  features  of  the  model  for  reflection  from  a  widespread 
disturbance  are  illustrated  in  Fig.  2,  showing  (from  the  top)  a  plane 
waveguide  mode  with  vertical  electric  field  E^  incident  at  angle  6  on 
a  disturbed  region,  which  is  assumed  stratified  in  one  lateral  dimen¬ 
sion.  Taking  the  direction  of  stratification  to  be  parallel  to  the 
x-axis  yields  no  loss  of  generality,  so  the  disturbance  has  lateral 
variation  only  in  the  y-directlon.  (Recall  that  vertical  variations 


Disturbed  region. 


2— Schematic  of  reflection  from  disturbed  region  with  diffuse, 
stratified  boundary.  Arrows  show  direction  of  propagation; 
electric  vectors  point  in  z-di recti  on. 
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in  the  ionosphere  are  accounted  for  in  the  waveguide  propagation 
parameter  S.)  Our  task  is  to  compute  the  reflected  signal  E  ,  which 
is  ignored  in  the  often-used  two-dimensional  WKB  treatment  of  ELF 
propagation. 

The  uniform  waveguide  propagation  constants  Sq  and  respectively 
characterize  the  undisturbed  and  disturbed  regions.  Lateral  gradients 
in  the  boundary  cause  the  partial  reflection  of  the  incident  signal 
from  the  disturbed  region.  Under  realistic  conditions,  disturbed 
regions  of  the  ionosphere  are  separated  from  undisturbed  regions  by 
boundary  regions  whose  widths  compare  with  a  reduced  ELF  wavelength  A/2n. 
Given  those  conditions,  we  cannot  use  either  of  the  two  usual  approx¬ 
imations — abrupt  boundary  or  slowly-varying  boundary — to  calculate  the 
reflection  coefficient.  Instead,  we  must  devise  a  model  that  fully 
accounts  for  the  finite — but  nonzero — gradients  In  the  diffuse  boundary, 
and  permits  an  analytic  solution  to  Eq .  (6).  The  formula 

(sj  -  Sq)  exp  [ (y  -  yn)/Ay] 

S  "  S0  +  -  1  +  exp  [  (y  -  yQ )  MyT~~  (13) 

has  the  desired  features:  namely,  S(y)  -*■  Sq  in  the  undisturbed  region; 
S(y)  -*  in  the  disturbed  region;  and  the  transition  from  Sq  to  S^ 
occurring  continuously  over  a  distance  Ay,  which  can  be  varied  to 
describe  moderate,  abrupt,  or  gradual  boundaries  centered  at  y  *  yQ. 

Before  deriving  the  lateral  reflection  coefficient,  we  consider 
the  types  of  disturbances  the  model  outlined  by  Fig.  2  and  Eq.  (13) 
can  represent.  The  assumption  of  a  stratified  boundary  applies  only 
if  the  boundary  curvature  greatly  exceeds  A/2ir.  Similarly,  the 
assumption  that  the  disturbance  is  of  semiinfinite  extent  in  the 
y-direction  ignores  reflections  from  the  rear  edge,  although  they 
must  in  principle  exist  on  a  sphere  (the  earth).  Given  that  ELF 
waves  are  thousands  of  kilometers  long  and  that  A/2ir  is  typically 
between  0.5  and  1  Mm,  our  model  thus  applies  only  to  disturbances 
that  cover  large  portions  of  the  earth  or  to  the  diurnal  effect  on 
the  waveguide,  where  the  diffuse  boundary  represents  the  terminator. 
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IV.  REFLECTION  COEFFICIENT 


This  section  derives  an  equation  for  the  reflection  coefficient, 
solving  Eq .  (6)  for  the  waveguide  model  in  Sec.  III.  Recall  from  Eq . 

(5)  that  E  [governed  bv  Eq .  (6)]  must  be  multiplied  by  the  local  excita¬ 
tion  factor  A  to  obtain  the  true  electric  field.  Both  the  incident  and 
reflected  waves  have  the  same  excitation  factor,  which  therefore  cancels, 
giving  for  the  reflection  coefficient 


E  E_ 

R-  J-  J-R 

1  EI 


However,  the  transmitted  wave  contains  A^  of  the  disturbed  region  of 
the  waveguide,  such  that  the  true  transmission  coefficient  is 


T  =  ^ 

E1  A0  Ej  A0 


The  difference  between  T  and  T  must  be  accounted  for  if  transmission 
coefficients  are  to  be  calculated  from  the  reflection  coefficients 
given  below.  Unlike  at  higher  frequencies,  higher  order  modes  generated 
by  the  boundary  need  not  be  considered  at  ELF,  where  all  but  the  TEM 
mode  are  far  below  the  cutoff  frequency. 

Substituting  Eq.  (13)  into  Eq .  (6),  and  considering  that 
E  <*■  exp  (-ikx  sin  9)  for  stratified  media,  we  obtain  the  following 
form  for  the  two-dimensional  wave  equation: 


A2  „  2  )  2  (N  -  1)  exp  [(y  -  yQ)My])  „ 

— ~  E  +  kn  i  cos  0  +  — — - P7 - r-7T — i - >  E 

d  2  01  1  +  exp  [ (y  -  y0)/AyJ  t 


J  E  -  0  . 


2  2  2 

For  Eq.  (16)  we  define  N  -  S^/Sq — physically  signifying  a  complex 

2  ?  o 

refractive  index— and  kQ  =  *  S“ — denoting  the  complex  wave  number 
in  the  undisturbed  portion  of  the  waveguide. 


I 
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To  calculate  the  reflection  coefficient,  we  apply  the  following 
changes  of  variable  to  convert  Eq.  (10)  to  a  standard  form: 

-ik  Aycos0£ 

E  =  e  (1  +  e5)MO  ,  (17) 

where  £  =  (y  -  y^)/Ay.  Considerable  but  straightforward  rearrangement 
yields  for  \J) 

(1  +  e^)tp"  +  [-2ik^  Ay  cos  0  +  (2  -  2ikg  Ay  cos  0)e^]i^' 

+  [1  -  2ikQ  Ay  cos  0  +  k2(Ay)2(N2  -  l)]eS  =  0  ,  (18) 

where  a  prime  denotes  differentiation  with  respect  to  £. 

The  form  of  Eq .  (18)  is  such  that  hypergeometric  functions  can  be 
used  to  derive  the  reflection  coefficient.  The  procedure  is  only  out¬ 
lined  here,  since  it  is  similar  to  that  originally  used  by  Epstein 
[1930] — and  more  recently  documented  by  Budden  [1961] — in  studying  the 
reflection  and  transmission  of  radio  waves  that  are  incident  on  the 
ionosphere  from  below.  The  derivation  follows  from  recognizing  that 
the  solutions  to  the  equation 

(1  +  eC)<J>"  +  [c  -  1  +  (a  +  b)eC)<J»'  +  abeSj)  =  0  (19) 

*  C  -1 

are  hypergeometric  functions  of  argument  n  ■  -e  and  n  that  con¬ 
tain  the  usual  coefficients;  that  is,  the  solutions  to  Eq.  (19)  are 
combinations  of  F(a,  b,  c;  n)»  F(a,  l-c+a;l-b+a;n  ),..., 
where  F  is  the  standard  hypergeometric  function  [Abramowltz  and  Stegun, 
1964]. 

Well-known  analytic  continuation  formulas  connect  the  hypergeo¬ 
metric  functions  of  argument  n  with  those  of  argument  r)  The 

By  the  change  of  variable  n  *  -e^,  Eq.  (19)  transforms  to 

n(l  -  n)  ^-9  +  [c(a  +  b  +  l)n ]  -  ab\(>  -  0  , 

dn  ari 

which  is  the  best-known  form  of  the  hypergeometric  equation. 


c 
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leading  terms  of  power-series  expansions  are  then  applied  in  both  the 
undisturbed  region,  where  n  <<:  1.  and  the  disturbed  region,  where 
r)  ^  «  1.  The  expansions  contain  readily  identifiable  terms  such  as 
E  ,  Er,  and  E^,  from  which  the  reflection  and  transmission  coefficients 
are  easily  formed.  The  steps  are  tedious  but  straightforward,  and 
lead  to  the  formula  for  the  reflection  coefficient  [Budden,  1961] : 


r(c  -  nrq  -  a)rq  +  b  -  c) 
r(c  -  a)r(b)ru  -  c) 


(20) 


where  T  is  the  gamma  function  [Abramowitz  and  Stegun,  1964]. 

The  final  step  relates  our  model  to  the  canonical  result  of  Eq. 
(20)  through  the  equivalence  of  Eqs.  (18)  and  (19) — which  obtains 
provided 


a  =  1  -  ikQ  Ay(cos  0  +  q^)  , 

b  *  1  -  ikQ  Ay(cos  0  -  q^)  ,  (21) 

c  =  1  -  2ikg  Ay  cos  8  , 

where 

qi  3  si/so  ~  sin2  6  *  (22) 

Substituting  Eq.  (21)  into  Eq.  (20),  and  using  the  identity  T(1  +  z)  * 
zT(z),  we  have 


where 


R  =  Ro(0)  D(0 ,  Ay)  , 
cos  0  -  q 

R„(6)  *  - q— .  -  1 

0  cos  0  +  q^ 


(23) 


(24) 


is  the  Fresnel  coefficient  for  reflection  from  a  sharp  boundary,  and 


-13- 


T(1  -  2ikSn  Ay  cos  0) [ 1  +  ikSn  Ay(q,  +  cos  9)] 

D(0,  Ay)  - - 5 - 2 - - - “ - (25) 

r(l  +  2ikSQ  Av  cos  9)F  [1  +  ikS^  Ay(q^  -  cos  9)) 

is  a  diffuseness  factor  that  accounts  for  the  reduced  reflection  due 
to  the  nonabruptness  of  the  boundary.  It  is  easily  shown  that  |d|  be¬ 
comes  smaller  as  Ay  increases.  Moreover,  when  Ay  =  0,  D  *  1,  and 
Eq.  (23)  reduces  to  the  Fresnel  coefficient  R^. 


t 


t 
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V.  NUMERICAL  RESULTS 


This  section  presents  numerical  results  for  the  modal  reflection 
from  a  disturbed  region  of  semiinfinite  extent  and  for  the  relative 
propagation  functions  for  laterally  bounded  disturbed  regions.  To 
emphasize  reflection  phenomena,  most  reflection  coefficients  are  cal¬ 
culated  for  45  Hz,  the  lowest  ELF  transmitting  frequency.  The  rela¬ 
tive  propagation  function  is  calculated  for  75  Hz. 

SEMIINFINITE  DISTURBED  REGION 

In  the  following  pages,  the  reflection  coefficient  R  [calculated 
from  Eq.  (23)]  is  plotted  as  a  function  of  incidence  angle  0  and 
boundary-thickness-parameter  Ay.  In  all  but  one  case  the  incident  and 
reflected  waves  are  assumed  to  be  in  the  ambient  daytime  region  of  the 
earth,  where  the  waveguide  propagation  constant  is  Sq  (see  Fig.  2). 

We  assume  SQ  =  1.2  -  0.09i  at  45  Hz  and  SQ  -  1.15  -  0.0851  at  75  Hz, 
corresponding  to  daytime  attenuation  rates  of  0.73  dB/Mm  and  1.15  dB/Mm. 
Those  figures  agree  reasonably  well  with  both  measured  and  calculated 
values  [e.g.,  Pappert  and  Mole r,  1974;  Ginsberg,  1974;  Galejs,  1972], 
although  SQ  depends  somewhat  on  propagation  direction,  latitude,  and 
season,  and  values  larger  or  smaller  than  those  here  have  been  reported. 

The  table  below  gives  the  values  of  used  to  calculate  the  re¬ 
flection  coefficient.  The  various  SPE  entries  are  averages  based  on 
numerous  calculations  performed  over  several  years  [e.g.,  Field,  1969]. 
(The  strong  SPE  model  also  applies  to  a  moderate  nuclear  disturbance.) 
The  value  for  the  nighttime  propagation  constant  is  rather  problematic, 
in  that  experimental  measurements  exhibit  much  scatter  [Bannister, 

1974;  White  and  Willim,  1974].  Moreover,  the  measured  nighttime 
attenuation  rates  tend  to  be  somewhat  larger  than  those  calculated 
from  the  accepted  model  ionospheric  height  profiles  [Pappert  and  Moler, 
1974].  Nevertheless,  both  measurements  and  calculations  show  that  the 
attenuation  rate  and  ReS  are  usually  smaller  at  night  than  in  the  day¬ 
time.  The  ambient  night  value  shown  for  in  the  table  is  therefore 
slightly  smaller  than  the  ambient  day  value  for  Sq  given  above. 
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VALUES  OF  S1 


Condition 

Frequency 

(Hz) 

Attenuation 
Rate,  a 
(dB/Mrn) 

S1 

Strong  SPE 

45 

2.4 

1.75-0.31 

75 

3.4 

1.50-0. 25i 

Moderate  SPE 

45 

1.6 

1.45-0.21 

Weak  SPE 

45 

1.3 

1.35-1.161 

Ambient  night 

45 

0.6 

1.15-0.0751 

Figures  3  through  6  plot  R  against  the  incidence  angle  of  an  ELF 
waveguide  mode  impinging  on  the  boundary  that  separates  the  ambient 
daytime  region  of  the  earth  from  each  SPE  in  the  table.  In  Fig.  7, 

R  is  plotted  for  a  mode  incident  on  the  day /night  terminator  from  the 
nighttime  side,  and  in  Fig.  8,  for  a  mode  incident  from  the  daytime 
side.  In  all  these  figures,  the  curve  RQ  is  the  Fresnel  reflection 
coefficient,  which  applies  if  the  two  regions  are  separated  by  a 
sharp  boundary  [see  Eq.  (23)].  The  remaining  curves,  which  account 
for  the  diffuseness  of  the  boundary,  result  from  multiplying  Rq  by 
the  diffuseness  factor  [Eq.  (25)].  Their  labels  quantify  the  effec¬ 
tive  boundary  thickness  Ay. 

The  reflection  coefficient  has  several  common  features  in  Figs. 

3  through  6  that  are  relatively  insensitive  to  SPE  strength.  First, 
it  is  small — less  than  a  few  tenths — for  all  but  the  most  oblique 
incidence  angles.  Second,  it  is  much  reduced  by  the  diffuseness  of 
the  boundary,  particularly  if  the  boundary-transition  width  Ay  is 
greater  than  X/2tt.  Third,  it  is  smaller  at  75  Hz  than  at  45  Hz,  al¬ 
though  its  sensitivity  to  the  transition  width  Ay  is  greater  at  the 
higher  frequency  (compare  Fig.  3  with  Fig.  4). 

Satellite  data  [Reagan  and  Watt,  1976]  show  that  the  energetic 
proton  flux  during  the  August  1972  SPE  was  nearly  constant  at  invariant 
latitudes  above  deg,  rolling  off  rapidly  outside  the  polar  cap  and 
|  diminishing  an  order  of  magnitude  over  3  to  5  deg  of  latitude.  Since 


r  r 
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ram  weak  SPE  (F  =  45  Hz) 
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1  deg  of  latitude  corresponds  to  approximately  100  km  of  lateral  dis¬ 
tance,  we  infer  that  Ay  300  to  500  km  reasonably  represents  an  actual 
SPE.  For  that  Ay,  significant  reflection  (|r|  >,0.7)  will  occur  if 
the  TEM  mode  is  incident  on  the  SPE  boundary  at  angles  greater  than 
80  to  85  deg  (see  Figs.  3  through  6). 

Such  oblique  incidence  angles  may  actually  occur  on  continental 
U.S. -North  Atlantic  ELF  paths  that  are  nearly  tangent  to  the  polar 
cap.  Just  outside  the  cap,  we  would  expect  interference  patterns 
between  the  direct  ELF  mode  and  that  obliquely  reflected  by  an  SPE. 
in  addition,  since  poor  transmission  occurs  for  incidence  angles  at 
which  the  reflection  coefficient  is  large,  receivers  on  paths  that 
obliquely  intersect  the  polar  cap  will  be  in  "shadow  zones"  during  an 
SPE,  with  signals  degraded  by  reflection  beyond  the  amount  due  to  the 
increased  a  associated  with  the  SPE  (see  the  table  presented  earlier). 

Our  discussion  considers  only  the  reflection  of  modes  incident 
on  the  SPE  from  an  undisturbed  region,  it  being  unlikely  that  an  ELF 
transmitter  will  ever  be  placed  within  the  polar  cap.  But,  since  a 
transmitter  must  operate  under  both  nighttime  and  daytime  conditions, 
we  consider  both  the  reflection  of  waves  incident  from  the  night  side 
of  the  terminator  and  those  from  the  day  side  (see  Figs.  7  and  8). 

From  data  on  diurnal  phase  shifts  for  north-south  VLF  paths  [Chilton, 
Crombie,  and  Jean,  1961],  the  earth-ionosphere  waveguide  makes  the 
transition  from  daytime  to  nighttime  propagation  conditions  in  60  to 
90  min,  during  which  time  the  earth's  surface  rotates  some  1500  to 
2500  km.  Hence,  the  thickness  of  the  day /night  terminator  is  greater 
than  the  several-hundred-kilometer  thickness  of  the  SPE  boundary. 

The  values  of  Ay  used  to  construct  Figs.  7  and  8  are  therefore  greater 
than  those  used  for  Figs.  3  through  6. 

Figure  7  shows  that  the  reflection  coefficient  for  night-to-day 
propagation  is  very  small  for  all  realistic  incidence  angles,  and  that 
the  terminator  would  only  slightly  affect  ELF  signal  strength.  Such 
is  as  expected,  because  the  difference  S^  -  S^  is  small  for  night-to¬ 
day  propagation  and  the  boundary  is  very  diffuse. 

The  results  for  day-to-night  propagation  (Fig.  8)  are  at  first 
surprising.  Although  R  is  small  (as  expected)  for  small-to- intermediate 


incidence  angles,  it  becomes  large  ( | R |  >0.7)  for  oblique  incidence — 
a  phenomenon  analogous  to  the  total  internal  reflection  that  occurs 
when  a  light  wave  propagates  from  glass  to  air.  For  the  assumed  day- 
to-night  propagation  parameters,  <  Sq  (see  p.  15),  and  the  reflec¬ 
tion  coefficient  becomes  large  when  the  critical  incidence  angle 

8^  Arc  sin  S^/Sq  «  75  deg 

is  exceeded.  In  fact,  only  the  complex  nature  of  Sq  and  prevents 

the  reflection  from  being  total  for  9  >  0  . 

c 

The  internal  reflection  shown  in  Fig.  8  will  usually  be  unimportant 
because  ELF  paths  tend  to  be  east-west  rather  than  north-south.  How¬ 
ever,  it  could  prevent  a  daytime  transmitter  in  the  northern  hemi¬ 
sphere  from  reaching  a  nighttime  receiver  in  the  southern  hemisphere 
if  the  incidence  angle  on  the  terminator  exceeded  the  critical  angle, 
depending  on  the  day  and  night  waveguide  parameters. 

BOUNDED  DISTURBANCE 

The  algorithm  used  to  solve  Eq.  (11)  is  the  subject  of  the  Appendix. 
The  numerical  results  given  below  illustrate  the  dependence  of  the  fields 
on  the  strength,  location,  and  extent  of  the  disturbance.  Although  the 
results  pertain  to  the  analytic  model  defined  in  Eq.  (12),  the  algorithm 
will  work  equally  well  for  any  spatial  dependence  S(x,  y),  provided 
only  that  S(x,  y)  -  Sq  effectively  vanishes  outside  a  rectangle  several 
megaraeters  on  a  side. 

All  results  pertain  to  a  frequency  of  75  Hz,  and  as  above,  Sq  = 

1.15  -  0.085i  characterizes  the  undisturbed  portion  of  the  waveguide. 

The  propagation  constant  at  the  center  of  the  disturbance  is  taken 
to  be 

Sj  =  A(l. 5  -  0.251  -  SQ)  +  SQ  , 

so  that  strong,  moderate,  and  weak  disturbances  can  be  modeled  by  let¬ 
ting  A  =  1,  2/3,  or  1/3.  For  A  ■  1,  assumes  the  value  given  earlier 
for  a  strong  SPE. 
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Figure  9  plots  W  versus  the  off-path  distance  y  for  a  strong 
(A  =  1),  longitudinally  confined  (Ax  =  0.5  Mm)  disturbance,  with 
curves  shown  for  several  effective  half-widths  Ay.*  In  Fig.  9a,  W 
is  computed  from  the  simple  nondeviative  WKB  formula 

(26) 

where  the  integration  path  is  the  straight  line  connecting  transmitter 
and  receiver.  Equation  (26)  is  shown  in  the  figure  to  give  the  in¬ 
tuitively  expected  result;  namely,  W  is  smallest  for  y  =  0,  where  the 
direct  path  intersects  the  peak,  of  the  disturbance,  and  approaches 
unity  for  y  >>  Ay--in  which  case  the  direct  path  effectively  misses 
the  disturbance. 

On  the  other  hand,  W  calculated  from  the  integral  equation  [Eq. 
(11)]  is  shown  in  Fig.  9b  to  exhibit  a  characteristic  pattern  on  the 
line  x  =  10  Mm:  a  maximum  for  y  =  0,  a  minimum  a  few  megameters  off- 
path,  and  increasing  with  off-path  distance.  That  counterintuitive — 
but  correct — full-wave  behavior  of  W  is  due  to  focusing  and  diffrac¬ 
tion — the  effects  of  which  are  ignored  in  the  nondeviative  WKB  formula. 

Focusing  occurs  because  wave  normals  tend  to  bend  toward  regions 
where  S(x,  y) — which  in  the  waveguide  formulation  is  analogous  to 
the  refractive  index — is  large.  The  disturbance  therefore  behaves 
as  a  converging  cylindrical  lens.  Diffraction  occurs  because  the 
effective  aperture  filled  by  the  lens  is  approximately  2Ay,  which  is 
smaller  than  or  comparable  to  a  Fresnel  zone  for  the  parameters  used 

* 

To  minimize  computation  time,  the  results  in  this  subsection 
are  calculated  for  a  longitudinal  half-width  Ax  =  0.5  Mm.  The  effect 
of  the  disturbance  is  therefore  modest,  and  W  is  within  20  to  30  per¬ 
cent  of  unity  for  all  cases.  A  larger  value  for  Ax  would  increase 
the  anomalous  attenuation  and  therefore  yield  a  W  much  less  than 
unity.  Except  as  noted,  the  choice  of  Ax  affects  the  results  only 
in  detail;  it  does  not  alter  the  overall  conclusions.  In  Fig.  9  and 
subsequently,  the  total  pathlength  for  y  -  0  is  10  Mm,  and  the  dis¬ 
turbance  is  centered  at  midpath  (x  *  5  Mm) . 


r  — — - TP 


WKB 


exp  (  -ik 
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[S(x\  y’)  - 


so]  dr 


,) 


(a)  WKB  approximation 


y  (Mm) 


(b)  Integral  equation 


Fig.  9~Ampl1tude  of  W  versus  off-path  distance  y  for  x  =  li 
x  *  5  Wn,  Ax  »  0.5  Wn,  and  strong  (A  *  1)  disturbance 
with  various  half-widths 
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here.  As  expected,  the  diffraction/focusing  pattern  becomes  much  less 
pronounced  as  2Ay  approaches  the  width  of  the  first  Fresnel  zone. 
Focusing  may  actually  overcome  the  anomalous  attenuation  through  the 
center  of  the  disturbance,  causing  |w|  to  exceed  unity  (see  Fig.  9b). 

Figure  10  is  analogous  to  Fig.  9,  except  that  the  curves  are 
parametric  in  disturbance  strength  rather  than  in  transverse  half¬ 
widths.  As  expected,  the  effects  of  lateral  focusing  are  most  pro¬ 
nounced  for  the  strongest  disturbance.  A  direct  path  that  intersects 
the  peak  of  the  disturbance  (receiver  at  y  =  0)  is  illustrated  in 
Fig.  11,  showing  the  dependence  of  W  on  the  transverse  half-width  Ay 
and  treating  disturbance  strength  parametrically.  For  reasons  dis¬ 
cussed  above,  signal  strength  W  tends  to  decrease  as  the  disturbance 
broadens  (Ay  increases)  and  strengthens  (A  increases). 

To  assess  the  relationship  between  W  and  the  transverse  dimension 
Ay  of  a  disturbance,  it  is  useful  to  define  a  transverse  shape  factor 
Q,  such  that 


Q  = 


(27) 


For  a  disturbance  centered  on  the  direct  path,  the  WKB  solution  ignores 

transverse  variations.  Thus,  Eq.  (27)  is  the  ratio  of  the  correctly 

computed  propagation  anomaly  to  that  computed  by  ignoring  transverse 
* 

gradients.  Since  the  WKB  solution  overestimates  the  anomaly  by 
accounting  for  only  the  peak  of  an  on-path  disturbance,  Q  is  every¬ 
where  smaller  than  unity. 

For  long  paths  and  very  weak  disturbances  satisfying  the  conditions 


k|SQ|x  »  1  , 
k2|s2  -  Sq| (Ax)2 


(28) 


(29) 


Strictly  speaking,  the  WKB  solutions  also  ignore  backward  reflec¬ 
tions.  However,  the  effect  on  the  received  signal  is  small  enough  to 
be  safely  ignored  here. 
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and 


«  1  . 


(30) 


Field  and  Joiner  [1979] — using  a  perturbation  method  equivalent  to 
the  first  Born  approxiroation--showed  that  W  has  the  limiting  form 


Q 


VS0 


it  kAy;  J 


■1/2 


(31) 


where  d,  the  maximum  half-width  of  the  first  Fresnel  zone,  is  given  by 


Equation  (31)  shows  that  Q  -*■  1  and  the  transverse  gradients  can  be 
safely  ignored  if  Ay  >  d,  so  that  the  peak  of  the  disturbance  effec¬ 
tively  fills  the  first  Fresnel  zone.  Conversely, 

|Q|  if  ^«1,  (33) 

giving  the  realistic  result  that  for  y  =  0,  the  ratio  of  the  actual 
attenuation  to  that  calculated  from  the  WKB  formula  (ignoring  the 
finite  transverse  gradients)  corresponds  to  the  fraction  of  the  first 
Fresnel  zone  filled  by  the  disturbance. 

Figure  12  shows  Q  computed  numerically  from  Eq.  (33)  for  three 
values  of  the  strength  parameter  A.  For  comparison,  the  limiting 
value  given  by  Eq .  (31)  is  also  plotted.  The  simple  limiting  form 
for  Q  is  remarkably  accurate  even  for  A  *  1  (a  very  strong  ionospheric 
disturbance) . 

To  illustrate  the  longitudinal  dependence  of  the  field.  Fig.  13 
plots  |w|  versus  x  for  a  strong  disturbance  of  longitudinal  half-width 

*Sf  Sq  implies  A  -►  0. 


t 


Fig.  13--Amplitude  of  W  versus  distance  along  direct  propagation 
path  for  strong  disturbance  (A  =  1 ,  y  =  0,  Ay  =  1  ttn) 
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Ax  =  0.5  Mm  and  transverse  half-width  Ay  =  1  Mm,  centered  on  the  direct 
path.  The  standing  wave  pattern  in  front  of  the  disturbance  is  caused 
by  reflections  from  the  severe  longitudinal  gradients;  it  is  much  less 
pronounced  for  larger  values  of  Ax.  |w|  gradually  increases  with  dis¬ 
tance  behind  the  disturbance,  inasmuch  as  focusing  causes  energy  to 
converge  on  the  plane  y  =  0.  For  certain  off-axis  positions  of  the 
receiver,  the  electric  field  slowly  decreases  with  distance. 


J 


APPENDIX 


Numerically  calculating  the  effect  of  a  localized  perturbation 
in  the  waveguide  reduces  to 

1.  Solving  Eq.  (11)  [p.  6]  in  the  disturbed  region. 

2.  Integrating  Eq.  (11)  over  the  disturbed  region  to  obtain 
W(x,  y). 

2  2 

Assuming  a  localized  perturbation,  S  -  S^  vanishes  outside  some  bounded 
region,  allowing  us  to  partition  the  region  into  a  grid  such  that  Eq. 
(11)  takes  the  form 


W(x,  y) 


dy'  K(x,  y,  x', 


y’)W(x\  y')  , 


(A.  1) 


where  denotes  n  squares  centered  at  x^,  y^  in  the  grid.  Figure 
A. 1  illustrates  the  partition  for  solution  within  a  disturbance  (shaded 
portion);  Fig.  A. 2,  the  partition  for  outside  a  disturbance. 


SOLUTION  WITHIN  DISTURBED  REGION 

Equation  (A.l)  may  be  rewritten 


W(x0,  yQ) 


=  1  + 


// 


(1=0) 


dx'  dy'  K(xq ,  yQ,  x' ,  y’)W(x',  y') 


+  y]  ff dy'  K(xo*  yo*  x’’  y,)w<x’’  y,) 


(A.  2) 


where  W(Xq,  Yq)  indicates  that  the  solution  point  resides  in  square  0. 
Assigning  every  square  the  index  0  in  sequence  makes  Eq.  (A. 2)  rig¬ 
orously  equivalent  to  Eq.  (A.l).  It  follows  from  Eq.  (A. 2)  that 
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n- 1 


1  +  ^  W(x1,  y t )  //*■  '  dy'  K(xq,  yf 


,  x\  y') 


W(x0,  yQ)  s- 


1  ~JJ dx'  dy'  K(xQ ,  yQ,  x' ,  y') 


(A. 3) 


Equation  (A. 3)  entails  the  approximation  that  for  a  sufficiently 

fine  grid,  W(x',  y')  «  W(x|,  y!),  the  value  at  the  center  of  any 

square,  and  may  hence  be  factored  out  of  that  integral.  The  kernel  K 

is  integrated  with  standard  quadrature  methods  over  squares  i  =  1  to 

n  -  1.  Square  0  contains  a  subgrid  to  accommodate  the  singularity  of 

(2) 

K  due  to  Hq  ’ (kS0r2) . 

The  algorithm  becomes  an  iterative  solution  of  Eq .  (A. 3),  where 
all  W  are  initially  set  to  one.  A  single  iteration  consists  of  n 
solutions  in  which  every  square  assumes  an  index  in  the  permutation 
cycle  0,  1,  2,  3,  ...  n  -  1.  To  evenly  distribute  error  throughout 
the  grid,  the  cycle  starts  at  different  places  in  different  iterations. 

Since  only  a  few  kernel  functions  need  be  evaluated,  the  algorithm 
is  easily  optimized  by  storing  kernels  rather  than  recomputing  them; 
the  kernel  factors  recur  both  from  iteration  to  iteration  (a  self- 
evident  recurrence)  and  within  each  iteration  (not  as  self-evident 
but  more  important). 

Arbitrary  convergence  criteria  may  be  used.  In  the  present  prob¬ 
lem,  iterations  stop  at  the  same  count  regardless  of  whether  we  use  a 
mean  square  or  a  maximum  inter-iteration  error  measure.  The  size  of 
the  grid  squares  must  be  carefully  chosen,  since  the  computation  time 
required  for  a  solution  varies  as  the  fourth  power  of  the  number  of 
squares . 

SOLUTION  OUTSIDE  DISTURBANCE 

The  second  part  of  the  calculation  is  straightforward.  With  a 
receiver  outside  the  disturbance,  there  is  no  problem  of  singularity 
l  of  the  kernel.  Equation  (A.l)  may  then  be  solved  with  simple  quadra¬ 

ture  techniques,  where  W  inside  the  disturbed  region  is  known  from  the 
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first  part  of  the  calculation.  Only  the  geometric  significance  of 
terms  in  the  kernel  changes;  distance  from  the  observation  point  to 
the  disturbance  is  now  the  most  important  factor,  whereas  in  the 
solution  within  the  disturbance,  the  squares  adjacent  to  square  0 
dominated. 
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